{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f6de0869",
   "metadata": {},
   "source": [
    "# Some aspects of performance in image processing\n",
    "\n",
    "Image processing on scientific data can take a long time, especially in settings where a series of images are acquired (time series for example). It is sometimes desirable to reduce the run-time of a given image processing pipeline, for example to make real-time image processing possible. \n",
    "\n",
    "In this tutorial, we cover several aspects of this specific task."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7484967f",
   "metadata": {},
   "source": [
    "## Measuring computing times and profiling\n",
    "\n",
    "If you want to optimize execution time, you first need to measure it. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe158174",
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage import data\n",
    "import numpy as np\n",
    "img = data.eagle()\n",
    "print(img.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82b187d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.imshow(img, cmap='gray')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d30ffd90",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit\n",
    "from skimage import filters\n",
    "img_median = filters.median(img, footprint=np.ones((7, 7), dtype=bool))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0df60508",
   "metadata": {},
   "source": [
    "The %%timeit magic has many options (see https://ipython.readthedocs.io/en/stable/interactive/magics.html#magic-timeit). It is also possible to use the ``timeit`` module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe7446da",
   "metadata": {},
   "outputs": [],
   "source": [
    "from time import time\n",
    "from skimage import filters\n",
    "t1 = time()\n",
    "img_median = filters.median(img, footprint=np.ones((7, 7), dtype=bool))\n",
    "t2 = time()\n",
    "print(t2 - t1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b861905",
   "metadata": {},
   "source": [
    "Now we know the time taken by a given function, but we don't know why it should take a long time to execute. We can use Python's line profile to dive into the inner code of a function, and to identify which lines of the function take longer to execute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e322361f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage import color\n",
    "img = color.rgb2gray((data.retina()[300:-300, 300:-300]))\n",
    "out = filters.hessian(img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb19534c",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2)\n",
    "ax[0].imshow(img)\n",
    "ax[1].imshow(out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50720b5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from line_profiler import LineProfiler, show_func\n",
    "profile = LineProfiler()\n",
    "profile.add_function(filters.hessian)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e997dd2",
   "metadata": {},
   "outputs": [],
   "source": [
    "out = profile.runcall(filters.hessian, img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8dddac4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "profile.print_stats()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "943b43e5",
   "metadata": {},
   "source": [
    "**Exercise**: what is the function taking the largest part of the execution time inside ``hessian``? Add this function to the line profiler and line-profile this other function while executing ``hessian``."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3324ec5b",
   "metadata": {},
   "source": [
    "It is of course possible to apply the line profiler to your own functions. This can be useful when defining for example an image processing pipeline with a succession of operations. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "221f09b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext line_profiler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2eee335",
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage import data\n",
    "from skimage.filters import threshold_otsu\n",
    "from skimage.segmentation import clear_border\n",
    "from skimage.measure import label, regionprops\n",
    "from skimage.morphology import closing, square\n",
    "from skimage.color import label2rgb\n",
    "\n",
    "def get_labels(img):\n",
    "    # apply threshold\n",
    "    thresh = threshold_otsu(image)\n",
    "    bw = closing(image > thresh, square(5))\n",
    "    # remove artifacts connected to image border\n",
    "    cleared = clear_border(bw)\n",
    "    # label image regions\n",
    "    label_image = label(cleared)\n",
    "    return label_image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c337fc5",
   "metadata": {},
   "outputs": [],
   "source": [
    "image = data.coins()\n",
    "\n",
    "labels = get_labels(img)\n",
    "plt.imshow(labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4eb6513d",
   "metadata": {},
   "outputs": [],
   "source": [
    "%lprun -f get_labels get_labels(img)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4e7fff7",
   "metadata": {},
   "source": [
    "scikit-image uses a lot of cython code, and it is also possible to profile cython code. Please see https://cython.readthedocs.io/en/latest/src/tutorial/profiling_tutorial.html for more information."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed7054ba",
   "metadata": {},
   "source": [
    "## Parallelizing execution by chunking the image array\n",
    "\n",
    "One way to speed up the execution is to use the different CPU cores. Most scikit-image functions use only one core, so the other cores are not working. If we have several images, we can use different cores for the different images (using for example ``joblib.Parallel`` or ``dask``). But we can also divide the image into different sub-images and execute the function on the different sub-images. This is done by ``skimage.util.apply_parallel``.\n",
    "\n",
    "**Warning**: most functions have boundary effects. Therefore you must use some overlap between the different sub-images in order to avoid artifacts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ef37db9",
   "metadata": {},
   "outputs": [],
   "source": [
    "img = data.eagle()\n",
    "print(img.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c88996df",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit -n 1 -r 1\n",
    "from time import time\n",
    "out = filters.median(img, np.ones((11, 11)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9e02317",
   "metadata": {},
   "outputs": [],
   "source": [
    "import joblib\n",
    "joblib.cpu_count()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84c4ee55",
   "metadata": {},
   "source": [
    "We need to install the optional dependency dask which is used by ``apply_parallel``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f5a0188",
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install dask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9b8920e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit -n 1 -r 1\n",
    "from skimage import util\n",
    "out = util.apply_parallel(filters.median, img, depth=7, extra_arguments=[np.ones((11, 11))])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fade82f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit -n 1 -r 1\n",
    "from skimage import util\n",
    "out = util.apply_parallel(filters.median, img, chunks=500, depth=7, extra_arguments=[np.ones((11, 11))])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3a92c51",
   "metadata": {},
   "source": [
    "**Exercise**: try different sizes of chunks and compare the execution times."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90db8bd5",
   "metadata": {},
   "source": [
    "## Changing the algorithm\n",
    "\n",
    "Often, the biggest performance gains are obtained by using a different algorithm. There is often a tradeoff between execution time and algorithm specificity, as exemplified by the example below about denoising. (Example adapted from https://scikit-image.org/docs/stable/auto_examples/filters/plot_denoise.html#sphx-glr-auto-examples-filters-plot-denoise-py)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c9883bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from skimage.restoration import (denoise_tv_chambolle, denoise_bilateral,\n",
    "                                 denoise_wavelet, denoise_nl_means, estimate_sigma)\n",
    "from skimage import data, img_as_float\n",
    "from skimage.util import random_noise\n",
    "from time import time\n",
    "\n",
    "original = img_as_float(data.chelsea()[100:250, 50:300])\n",
    "\n",
    "sigma = 0.155\n",
    "noisy = random_noise(original, var=sigma**2)\n",
    "\n",
    "t0 = time()\n",
    "img_tv = denoise_tv_chambolle(noisy, weight=0.1, channel_axis=-1)\n",
    "t1 = time()\n",
    "print(f\"TV: {t1 - t0} s\" )\n",
    "\n",
    "t0 = time()\n",
    "img_nlmeans = denoise_nl_means(noisy,\n",
    "                channel_axis=-1)\n",
    "t1 = time()\n",
    "print(f\"Non-local means: {t1 - t0} s\")\n",
    "\n",
    "t0 = time()\n",
    "img_bil = denoise_bilateral(noisy, sigma_color=0.05, sigma_spatial=11,\n",
    "                channel_axis=-1)\n",
    "t1 = time()\n",
    "print(f\"Bilateral: {t1 - t0} s\")\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(nrows=1, ncols=4, figsize=(8, 5),\n",
    "                       sharex=True, sharey=True)\n",
    "\n",
    "plt.gray()\n",
    "\n",
    "# Estimate the average noise standard deviation across color channels.\n",
    "sigma_est = estimate_sigma(noisy, channel_axis=-1, average_sigmas=True)\n",
    "# Due to clipping in random_noise, the estimate will be a bit smaller than the\n",
    "# specified sigma.\n",
    "print(f'Estimated Gaussian noise standard deviation = {sigma_est}')\n",
    "\n",
    "ax[0].imshow(noisy)\n",
    "ax[0].axis('off')\n",
    "ax[0].set_title('Noisy')\n",
    "ax[1].imshow(img_tv)\n",
    "ax[1].axis('off')\n",
    "ax[1].set_title('TV')\n",
    "ax[2].imshow(img_nlmeans)\n",
    "ax[2].axis('off')\n",
    "ax[2].set_title('Non-local means')\n",
    "ax[3].imshow(img_bil)\n",
    "ax[3].axis('off')\n",
    "ax[3].set_title('Bilateral')\n",
    "\n",
    "\n",
    "fig.tight_layout()\n",
    "\n",
    "plt.show()\n"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "notebook_metadata_filter": "all,-language_info",
   "split_at_heading": true
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "orphan": true
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
